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Abstract 

Turing models have been proposed to explain the emergence of digits during limb development. However, 
so far the molecular components that would give rise to Turing patterns are elusive. We have recently 
shown that a particular type of receptor-ligand interaction can give rise to Schnakenberg-type Turing 
patterns, which reproduce patterning during lung and kidney branching morphogenesis. Recent knock- 
out experiments have identified Smad4 as a key protein in digit patterning. We show here that the 
BMP-receptor interaction meets the conditions for a Schnakenberg-type Turing pattern, and that the 
resulting model reproduces available wildtype and mutant data on the expression patterns of BMP, its 
receptor, and Fgfs in the apical ectodermal ridge (AER) when solved on a realistic 2D domain that 
we extracted from limb bud images of Ell. 5 mouse embryos. We propose that rcceptor-ligand-based 
mechanisms serve as a molecular basis for the emergence of Turing patterns in many developing tissues. 
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Introduction 

The limb has long served as a paradigm to study organ development. Its easy access enabled devel- 
opmental biologists to use tissue grafting to define organizing centers, and advanced genetics has led 
to the identification of the key regulatory circuits [I , ]. Early grafting experiments showed that tissue 
from the posterior limb bud (referred to as zone of polarizing activity, ZPA) induced ectopic digits when 
implanted in the anterior limb bud [3]. These extra digits formed a mirror image of the normal forelimb, 
i.e. 1234554321 in mice and 234432 in chicken. Here numbers refer to digits with the lowest number 
referring to the most anterior digit in the wildtype (i.e. the thumb in humans). The effect was shown 
to be concentration-dependent in that a smaller number of grafted cells induced fewer digits and failed 
to induce the most posterior identities in the anterior limb bud [ ]. Wolpert's French Flag model ex- 
plained these experiments by suggesting that digit number and identity were determined by the local 
concentration of a chemical compound [4]. This compound was predicted to be produced in the ZPA 
and to diffuse across the limb bud, thereby creating a concentration gradient. Decades after the grafting- 
experiments had suggested that patterning was induced by a diffusible compound, Sonic hedgehog (SHH) 
was isolated and shown to affect patterning in a concentration-dependent manner [5,6]. In agreement 
with the morphogen model, SHH is expressed in cells of the ZPA and diffuses across the limb bud [ ]. 
In mice that do not express SHH only digit 1 develops [8], and in the polydactylous mouse mutant extra 
toes (Xt) SHH is expressed also in the anterior region of the limb bud [9]. 

However, several other experiments question this simple morphogen model and suggest that a more 
complex mechanism determines digit patterning during limb bud development. For one, Polydactyly is 
still observed even when SHH expression is abolished in the extra toe (Xt) mutant [10]. Moreover, an 
experiment in which SHH expression is stopped prematurely at different stages during limb development 
results in an alternating anterior-posterior sequence of digit formation that cannot be reconciled with 
the standard morphogen model [ ]. Thus the shortest pulse that yielded an additional digit induced 
formation of digit 4, followed by digit 2, and then digits 3 and 5. Based on the morphogen model digits 
would be expected to emerge simultaneously or in an anterior-posterior sequence. 

Experiments in which SHH-expressing cells were marked revealed that descendants of the SHH- 
producing cells could be found far beyond the small posterior zone in which SHH expression can be 
detected [7, 12]. Descendants of the SHH-producing cells filled almost half of the limb bud and encom- 
passed all cells in the two most posterior digits and also contributed to the middle digit [7, 12]. Based 
on this observation it was suggested that it is the length for which cells express and secrete SHH (and 
therefore experience a high SHH concentration) that determines digit specification [ , ]. This model, 
however, still cannot explain the observed sequence of digit formation. 

A number of mathematical models have been developed to explain limb bud development and digit 
patterning. Based on experimental observations, the earliest model defined some 40 years ago a set of 
simple rules to recapitulate the observed changes in the limb buds' size and shape [13]. The modelling 
effort highlighted the importance of oriented cell division to obtain the elongated shape that is charac- 
teristic of the limb bud. Recent 3D imaging data in combination with a Navier-stokes description of the 
tissue as an viscous fluid indeed confirmed that isotropic proliferation alone is insufficient to explain the 
shape and size changes during limb bud development [14]. Local gain rates much higher than the mea- 
sured proliferation rates as well as loss rates had to be included; these may reflect oriented cell division 
and directed cell migration. The description of the limb tissue as a viscous fluid was first introduced 
by Dillon and Othmer who studied a model where the local growth rate (the source term in the Navier 
stokes equation) depended on the concentration of two signaling factors, SHH and FGF [ ]. A more 
recent model for limb bud growth explores the impact of differential tissue elasticity [16]. Other models 
have focused on particular regulatory interactions that pattern the limb expression domains, i.e. the 
AER-ZPA regulatory interaction [17], the establishment of the proximal and distal signaling centers [18], 
and the diffusivity of SHH [19]. 

Most models, however, have focused on the skeletal patterns that emerge during limb bud development 
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[20-29]. A number of different model types have been studied ranging from continuous to cell-based 
models and from rule-based stochastic models to reaction-diffusion models; for reviews see [30,31]. Turing 
mechanisms are a particular type of reaction-diffusion models [32] and have been frequently invoked to 
explain the emergence of pattern during development [33]; they have been studied in the context of 
limb patterning since 1979 [34]. In Turing models the number of modes (patterns) increases with larger 
domain size. Newman and co-workers showed that the different number of skeletal elements in the 
stylopod, zeugopod, and autopod in anterior-posterior direction can in principle be achieved as a result 
of the increase in the proximal-distal length; here they have to assume that patterns emerge sequential 
in stylopod, zeugopod, and autopod and are fixed (frozen) upon their emergence [23]. Newman and 
co-workers further showed that parameter values can be identified that give realistic steady state pattern 
if the equations are solved on static domains whose shapes are based on the limb shapes of various 
fossils [35]. However, there is no evidence that the patterning domain in the embryo has the same shape 
as the adult fossil and patterning is well known to occur on a growing rather than static domain. Miura 
and co-workers explored how the speed of pattern emergence could be affected [24] and studied Turing 
patterning on growing limb domains [22]. In particular, Miura and co-workers show that a Turing model 
can explain the increased number of digits in the doublefoot mutant with the observed increased domain 
size [22]. 

The many computational studies show that Turing models, in principle, are sufficiently flexible to 
account for the experimentally observed pattern in wildtype and mutant limb buds. However, in spite of 
the demonstrated sensitivity of Turing models to domain size and particular reaction types and kinetics 
the models have so far not been solved on realistic domains with biochemically realistic reaction kinetics. 
Recent efforts in this direction have been made by the Newman group who propose TGF-beta family 
ligands as the activator and FGFs, Noggin, Notch and/or CHL2 as the inhibitor in their Turing models 
with Schnakenberg kinetics [21,23,35]. However, mutant mice with defects in TGF-beta [36], FGF [3 —39], 
Noggin [40], and Notch [41] signaling all have digits, while CHL2 expression has been found to be restricted 
to chondrocytes of various developing joint cartilage surfaces and connective tissues in reproductive 
organs [ ]. While it is possible that the apparent robustness of the patterning mechanism to these 
mutations is a consequence of redundant regulatory interactions, we wondered whether the reported 
regulatory interactions would allow us to construct a model that would be consistent with available 
experimental data and that would not have to rest on missing data and redundancy arguments. 

Advances in experimental techniques now provide us with detailed information regarding the domain 
geometries [14], the timing of the processes, and the regulatory interactions [1,2,43]. Limb buds first 
form around embryonic day (E)9.5 while expression of Sox9, an early marker for tissue condensations, 
is first observed around E10.5 in the central mesenchyme of the limb bud (Fig. 1A, [44]). Around Ell 
the Sox9 expression pattern splits and a region in the stalk patterns the stylopod, two patches mark the 
zeugopod, and a region of strong expression marks the autopod (Fig. 1B,F, [44,45]). Around Ell. 5 the 
autopod expression pattern splits into distinct patches that mark the different digits (Fig. 1C,G, [44,45]). 
Over the next day more defined and elongated Sox9 expression pattern emerge that resemble the final 
digit structure (Fig. 1D,E,H [44,45]). Sox9 expression has been shown to be stimulated by BMP-2 [46] 
and the BMP-2 dependent increase in Sox9 then triggers increased expression of Noggin [4 7] which acts 
as an antagonist of BMP signaling by sequestering BMP in inactive complexes [48]. Noggin null mice 
form digits, but the digit condensations are much wider than in wildtype [40]. Much as Sox9, Noggin is 
expressed in digit condensations and serves as a marker of endochondral differentiation [40]. 

Even though BMPs also induce the expression of the BMP antagonist Gremlin [49] , Gremlin (unlike 
Noggin) is expressed in the interdigital space rather than in the digit condensations [50]. Another study 
notes a concentration-dependent regulation of Gremlin expression in that low concentrations of BMP2 
upregulate Gremlin while high concentrations of BMP2 downregulate Gremlin in limb mesenchyme cul- 
tures [51]. It is thus possible that stronger BMP signaling in the digit condensations induces Sox9 while 
weaker BMP signaling in the interdigital space induces Gremlin. Similarly, it has been suggested that 
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low BMP concentrations support expression of FGF8 while high BMP concentrations repress FGF4 and 
FGF8 expression in the apical ectodermal ridge (AER) [52,53]. Two different BMP receptors are ex- 
pressed in the limb mesenchyme, BMPRla (ALK3) and BMPRfb (ALK6) that bind the different BMP 
ligands with different affinities. BMPRla binds mainly BMP4 while BMPRfb binds BMP4, BMP2 and 
less well to BMP 7 [ , ] . The BMP4-BMPRfb affinity has been established in a Scatchard plot as Kd = 
530 pM [55]. Activin receptors (ActRI/ALK2 and ActRII) are also expressed in the limb mesenchyme [56] 
but appear to bind their BMP ligands with lower affinity [57]. BMP7 binds ActRII with 3-fold higher 
affinity compared to BMP2 (Kd = f .7 fjM versus 5.4 /LtM) and binds ActRI/ALK2 only with low affinity 
(Kd = 143 /LtM) [ ] . ft should be noted that the values are based on Biacore measurements and thus do 
not adequately reflect the cooperative binding effects. 

In both mouse and chicken expression of BMPRfb is restricted to precartilaginous condensations 
[58,59] while Bmprfa is expressed at low levels throughout the limb bud mesenchyme [59,60]. In the 
chicken Bmprla expression is elevated at the border between precartilaginous condensations and mes- 
enchyme, but low or absent in the precartilaginous condensations. Accordingly, overexpression of dom- 
inant negative (DN)-BRK-2 (the chicken receptor that corresponds to the mouse BMPRfB receptor) 
blocks chondrogenesis, but not overexpression of (DN)-BRK-f (which corresponds to mouse BMPRf A) 
or DN-ActR, 1 [61]. BMPR1B thus appears to be the major transducer of BMP signals in chicken limb 
condensations. 

In the chicken the two receptors have been shown to be engaged in different regulatory interactions. 
Thus BMP beads induce noggin, tgft>2, and bmprlb expression only when implanted into the tip of the 
chicken digit, but not when implanted in the interdigital space. Also bmprfa is not expressed when BMP 
beads are implanted in the interdigital space while expression of the bmpR-lb gene was enhanced in the 
adjacent digits with the spatial distribution appearing displaced toward the margin of the digits adjacent 
to the bead [62]. It therefore seems that in the chicken only BMPRlb (but not BMPRf a) enhances its 
own expression in response to ligand binding. TGF-beta beads also induce the expression of noggin in 
the interdigital regions but only with a time delay: noggin expression is detected after 30 hours, while 
bmprfb expression becomes detectable in the interdigital region already after 10 hours. TGF-beta thus 
most likely induces noggin expression indirectly by inducing bmprlb expression [62]. 

There are important differences between mouse and chicken limb development as there are differences 
between forelimb and hindlimb development. We will therefore mainly focus on mouse forelimb bud 
development and comment on other limb buds as appropriate. In the mouse both Bmprla and BmpRlb 
expression is elevated in the precartilaginous condensations and Bmprla and Bmprlb can compensate 
for each other [59,60]. BMPRlb null mice show normal digit patterning and no differences in Sox9 
expression are observed in wild-type and BmprIB-/- limbs up to and including Ef 2.5, indicating that the 
prechondrogenic limb mesenchyme is specified and is able to form condensations up to this time [60,63,64]. 
BmprlA-/- embryos are not viable and the severity of the BmprIA-/- conditional knock-out depends on 
the stage at which BmprIA-/- is removed. Cre-mediated excision in Col2-Cre mice seems to occur only 
in forming cartilaginous condensations, while in Prxf-Cre mice, it acts in limb mesenchyme at much 
earlier stages (95% active by E.10.5) [63]. If BMP receptor type IA is removed with Col2-Cre digit 
condensations can still be observed while in Prxl-Cre mice no digits form in the forelimb [ 13]. The 
mutants have shortened limbs and show almost complete agenesis of the autopod because of reduced cell 
proliferation [63]. Moreover, the expression of BMP target genes Msxf/2 and Grem are severely reduced. 
BMPRfb appears to rescue Msxl/2 expression. However, the expression is then colocalized with the 
patchy BMPRIB expression [63]. In the mouse BMPRla thus appears to be the major transducer of 
BMP signals in limb condensations. 

In the limb there are three important BMP ligands: BMP2, BMP4, and BMP7 [65]. BMP2 and BMP4 
appear to be the most important BMP ligands, but also in the BMP2: BMP4 double conditional knock- 
out some digits still form [66]. At E12.5 BMP2 is expressed predominantly in the interdigital webbing, 
BMP4 is expressed along the whole AER, and BMP7 expression is ubiquitous in the limb bud [65]. 



■5 



Bmp2 conditional mutants form five digits, conditional inactivation of a gene for BMP4 results in a 
polydactylous phenotype [67], but in BMP2:BMP4 conditional mutants two posterior digits are missing 
even though the limb field is broader [66]. Bmp7 null mice show an occasional anterior Polydactyly [65]. 
Importantly, even though BMP4 expression is reduced in the conditional BMP4 knock-out BMP signaling 
activity is increased [68] . To avoid the redundant role of multiple BMP ligands and receptors Zeller and 
co-workers conditionally removed the downstream signaling protein Smad4 (Co-Smad). The mutants still 
express Sox9 in the autopod but the pattern does not split into digits [69]. Smad4-dependent signaling 
is therefore necessary for digit formation. 

BMP signaling is embedded in a larger network that most importantly comprises FGFs and SHH in the 
limb. BMP signaling interferes with FGF-dependent signaling, at least in part by repressing expression 
of its receptor FGFR1 [ ] . FGF- loaded beads in turn repress BMPRlb and Noggin expression in chicken 
limb buds [ ]. In mouse limb buds FGFs repress the expression of the BMP antagonists Gremlin [71] and 
enhance the expression of SHH in the ZPA, which in turn enhances Gremlin expression. SHH induces the 
expression of BMP-2 and BMP-7 in chicken limb buds, and with a delay that of BRK2 (BMPRlb) [61]. 
Without SHH digits cannot form unless expression of Gli3 is removed as well [10]. In the absence of SHH 
signaling Gli3 forms the Gli3 repressor which prevents the expression of Gremlin and many other genes. 
However, in the absence of both SHH and Gli3 more digits appear than in the wildtype. Neither SHH 
nor Gli3 are therefore necessary for digit formation. Both Fgf4 and SHH expression terminate around 
E12 [43] but Fgf8 continues to be expressed in the AER. 

Given the complex regulatory interactions it is difficult to understand and predict the regulatory 
outcome by verbal reasoning. We therefore sought to build and analyse a computational model of these 
interactions. In light of the presence of digits without Shh/Gli3, and the absence of digits in BMPR1A and 
SMAD4 conditional mutants we focus on BMP-dependent signaling to explain the observed patterning 
of Sox9 in the limb bud. We find that the interactions between BMP and BMP receptor give rise to a 
Schnakenberg-type Turing mechanism. When we solve the model on a realistic domain we obtain patterns 
that are similar to those observed in experiments with wildtype and mutant mice. We propose that the 
observed patterning during limb bud development may be the consequence of a Turing type patterning 
mechanism that arises from the interactions between BMPs and their receptor. 

Results 
Model 

The core signaling pathways that regulate limb patterning have been well characterized and comprise, 
SHH, GLI, Gremlin, BMP, FGF, and their receptors [1,72]. Digit patterning is still observed in the 
absence of Gremlin [73], Noggin [40], SHH and GLI3 [10], but not in the absence of BMPR1A [63] and 
Smad4, an essential protein in BMP signaling [69]. We therefore neglect all dispensable network elements 
and focus on BMP signaling to explain digit patterning. We will include FGF as an important modulator 
and growth factor. Shh-dependent signaling is also an important modulator of BMP signaling in the limb 
bud and is clearly important to establish the anterior-posterior polarity. Turing patterns have been shown 
to persist in the presence of such a modulating gradient [74]. Accordingly, we focus on the mechanism 
that enables the emergence of digits in the autopod and we will neglect the detailed effects of Shh and 
how these enable the establishment of anterior-posterior polarity that results in the formation of a thumb 
and a pinky. 

BMP (which we denote by B) and FGF (denoted by F) are secreted proteins and diffuse much faster 
than the BMP receptor (denoted by R) which resides in the plasma membrane. Experiments on the 
BMP Drosophila homologue Dpp suggest that the majority of ligand-bound receptors, C, are internalised 
rapidly and reside mainly inside the cell [75]. We will therefore ignore diffusion of the receptor- ligand 
complex. We write Dj (j = B,R,F) for the diffusion coefficients with <C D^,Dp [76,77]. We write 
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D.A[-] for the diffusion fluxes where A denotes the Laplacian operator in Cartesian coordinates, and [•] 
concentration. The characteristic length of gradients depends both on the speed of diffusion and their 
removal. In the absence of contrary experimental evidence we will assume the simplest relation, linear 
decay, at rates dj[j] for all components (i.e. j = B, R, C, F). 

BMPs are dimers and one BMP molecule can therefore bind two receptors [78]. The rate of BMP- 
receptor binding is therefore proportional to R 2 B, and we use k on and k // as the binding and dissociation 
rate constants. BMP2 signaling has been shown to reduce BMP2 expression in the limb bud [79] and 
we therefore make the rate of BMP production negatively dependent on the BMP-receptor complex, C, 
i.e. we write for the BMP production rate Pb Kb+[c] anc ^ *^ us obtain f° r t ne BMP and BMP-receptor 
dynamics 

[B] = £BA[Bj+-PB ^-^p ^b[B] -k on [R] 2 [B]+fc off [C] (I) 

diffusion "v y degradation complex formation 

production 

[C] = k on [R] 2 [B]-fc off [C] -d c [C] . (2) 

complex formation degradation 

Signaling of BMP-bound receptors positively regulates receptor production [62], and we therefore must 
make receptor production dependent on the concentration of C. In the absence of contrary data we will 
use the simplest possible relation for the receptor production rate and write Pr([C]) = pr + pc[C] where 
Pn and pc are constants, 

[R] = D R A[R}+ p R +pc([C}) -d R [R] -2k on [R] 2 [B] + fc off [C] . (3) 

diffusion production degradation complex formation 

We can simplify this set of equations if we assume that the dynamics of the receptor-ligand complex are 
fast such that the receptor-ligand complex assumes a quasi steady-state. The concentration of bound 
receptor, C, is then proportional to R 2 B, i.e. 

[C] ~ r fc °" [R} 2 [B}=K c [Rf[B}; K c = - Kn . (4) 

k ff + do k ff + d c 

Equations (1,3) are sufficient for pattern formation and would reduce to the classical Turing model of 
Schnakenberg-type if we were to set pc = 2dc, and ds = 0. 

BMP expression and activity is modulated by FGF and Gremlin [68]. BMP induces Gremlin ex- 
pression and Gremlin then binds and sequesters BMP in an inactive complex [51,80]. Gremlin can in 
principle be secreted, but experimental evidence suggests that BMP4 activation and secretion are also 
negatively regulated by an intracellular Gremlin-BMP4 interaction [81]. We therefore do not consider 
Gremlin diffusion. Since experiments suggest that also Gremlin is a dimer [82] we assume 1:1 binding 
between BMP and Gremlin at rate fco n [B][G] and dissociation at rate k' Q ^[BG\ where we write G for 
Gremlin and BG for the BMP-Gremlin complex. If we model BMP-induced expression by a Hill function 
w ith Hill constant Kq and Hill coefficient n, and assume that dissociation of the complex (at 
rate k' Q ^) is much faster than its degradation then the previous system of equations would need to be 
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expanded to include 

Ami -i- P r . 

' K B + [C] 



I->] : MiDj-- ~dy[B] -k; n [R] 2 [B] + fc; ff [C] -( k; n [G][B]-fc; ff [BG] ) 

degradation complex formation complex formation 




production 

m = p, ln] ^l r k °"[ G ][ B ] + fc off[ BG l 

~* degradation complex formation 

production 

[BG] = k; n [G][B]-^ ff [BG]. (5) 

v v ' 

complex formation 

Experiments show that Gremlin is induced by BMP already early during limb development [68]. If we 
assume rapid formation of the BMP-Gremlin complex in comparison with other reactions we can use a 
quasi-steady state approximation for complex concentration: 

[BG] = kUG][B]-fco ff [BG]=0. (6) 

v ' 

complex formation 

The terms for complex formation thus disappear from the equations for Gremlin and BMP and the 
equation for Gremlin uncouples from the other equations. Gremlin seems thus not to affect the patterning 
module directly. We therefore ignore the Gremlin-BMP interaction in this model, and we will show that 
the observed oligodactily in the Gremlin knock-outs can be explained with the reduced limb bud size. 

Given the importance of FGFs for limb patterning we extended the model to also include FGF 
signaling. Much as BMP, FGFs can diffuse rapidly [83]. BMP expression is induced by FGF and the 
parameter P B in (1) thus becomes a function of FGF, F, i.e. Pb(F) = Pb + P*b [fv^+k 71 Kb+[C] ' wnere 
Kbf and n and the Hill constant and Hill coefficient respectively. Signaling of BMP-bound receptors, 
C ~ R 2 B (4), has been found to both stimulate and inhibit FGF-dependent processes [53,84]. Such 
inconsistent observations may be the result of a concentration-dependency of signaling, and it has been 
suggested that high levels of BMP signaling interfere with FGF-dependent processes while lower levels 
may positively regulate FGF expression in the AER [ ]. Hence FGF activity is best described as 
Pf([C]) — Pf jc^+k" [C]"+K n ' wri ere Kp\ -C Kp2 are the Hill constants for the activating and the 
inhibitory impacts of BMP signaling and n is the Hill coefficient. Here we note that FGF expression is 
restricted to the apical ectodermal ridge (AER) , and pf is therefore non-zero only on the distal boundary 
of the domain (dashed part of Fig. 2C; simulated expression patterns are shown in Fig. SI, S3G-J). We 
then write for the FGF dynamics 

[F] = D E ^ + „ w ll Jtj . w ^- jjjjj . (7) 

diffusion ' v / degradation 

production 

FGFs repress the expression of the BMP antagonist Gremlin, and enhance SHH expression in the ZPA, 
which in turn enhances Gremlin expression. However, as discussed above the Gremlin dynamics uncouple 
from the other equations, while the expression of SHH terminates before the expression pattern of Sox9 
has assumed the characteristic digit pattern [43-45] and digit patterning can still be observed in limb 
buds when expression of Shh (and Gli3) is removed in knock-outs [10]. We thus neither include Gremlin 
nor SHH in the model. 
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Summary We propose that the observed patterning in BMP signaling is the result of regulatory in- 
teractions between BMP, B, its receptor R, and FGFs, F as summarized graphically in Fig. 2A. The 
effective interactions in this model are graphically summarized in Fig. 2B. As final set of equations we 
have 

[Bl = D B A[B1 +p b +p* R r n 1 J B r 1or , ~d c Kc\R} 2 \B} - d B [B] 

ii 1 Fb PB [F] n + K" F K B + K C [R] 2 [B] m J 

[7?] = D R A[R]+p R + (pc-2dc)K c [R] 2 [B} - d R R 

rj* 75 A[pl , n (Kc[R] 2 [B]r K% 

~ DFA[F]+PF (K^mW+Kfi (K c [R?[B}) n + K- F2 ~ dF[n (8) 

We note that FGFs are not part of the core patterning mechanism. For large FGF concentrations 
([F] S> Kbf) the equation for F uncouples from the other two equations and if the constitutive BMP ex- 
pression pb, p*b was sufficiently large patterning could be obtained also without FGFs. Further regulatory 
interactions (like those with Gremlin) may tune the parameter values. 

Parameter Values The measurement of the parameter values (i.e. diffusion coefficients, production 
and degradation rates) in vivo is complicated and has only been carried out in few model systems, but 
not in the limb [86-88]. However, our conclusions do not depend on the exact values of parameters, 
but mainly depend on their relative values as can be seen by non-dimensionalizing the model. To non- 
dimensionalize the model we need to choose characteristic length and time scales, as well as characteristic 
concentrations. As characteristic time scale we use T = 1/Sf, i.e. t — Tt. As characteristic length 
scale we chose the characteristic length of the FGF gradient A = ^Dp/dp, i.e. x — XC,. We non- 
dimcnsionalize the FGF concentration with respect to the Hill constant Kbf, i-e. F — [F]/Kbf, and 

1 /3 

the BMP and receptor concentrations with respect to the characteristic concentration c = (K R /Kc) , 
i.e. B = [B]/c , R = [R]/cq. After re-grouping parameters as Si = Jj- except 5c = ^K c Cq, pi = 

except pp = v = %K G c\, Ki =^ = ^-, and D B = D R = §f we obtain 

3B F n 1 

— = D B A C B + PB+ pB l + Fnl + R2B - 5 C R 2 B S B B (9) 



D R A C R + p R + (p - 2S C )R 2 B ~ S R R (10) 



dR 

dF (R 2 B) n k™ 

a7 ~ c + pF (R 2 B) n + jq (R 2 B) n + ~tq ~ ' 

The dimensionlcss model contains five parameters less and the patterning mechanism no longer depends on 
absolute diffusion and decay constants, but only on the relative diffusion coefficients and the relative decay 
rates. Similarly, the absolute protein concentrations do not matter but only the relative concentrations 
(as result from the relative expression and decay rates) relative to the Hill constant and the effective 
binding constant Kq- 

Domain and Boundary conditions The equations were solved on a domain whose shape was ex- 
tracted from limb bud images at E12.5 as well as on a growing domain. The idealized domain is shown 
in Fig. 2C. We used zero-flux boundary conditions for B and R as the developing limb does not ex- 
change with the surrounding except at the flank which is ignored here. FGF production in the AER was 
implemented as a flux boundary condition, i.e. 

{R 2 B) n tq 

n - WF - P F (R2 B )n + K n {R 2 B) n +K n ^ 
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where ft refers to the unit normal vector. We also considered alternative implementations with production 
in a thin layer (Fig. S1B-E) or with Eq. 11 defined on the boundary (Fig. S1F). A thin layer with 
flux boundary conditions (Fig. SIC) and with Eq. 11 defined on the boundary (Fig. SIC) give similar 
patterns to those observed with flux boundary condition on the outer boundary (Fig. SlA). A thin layer 
with either constant (Fig. S1E) or BMP-dependent FGF production (Fig. SID) yielded overall the same 
pattern (same number of spots) but the shape of the spots was different. 



Modelling Domain Growth Proliferation is not uniform within the limb domain but the velocity 
field has not been published [14]. We therefore measured the proximal-distal axis of limb buds over 
developmental time and find that the growth rate is approximately linear (Fig. 3A) and we are therefore 
using a linear growth law with growth rate v g . We can then write for the proximal-distal length of the 
limb bud L(t) = Lq + v g t where Lq denotes the initial proximal-distal length of the limb bud, and t 
denotes time. 

To conserve mass (rather than concentration) the reaction-diffusion equations for the soluble factors 
B and F, equations 9 and 11 must be expanded to include the advection and dilution terms [89], i.e. 

® + V(n[X]) = D X A[X] + R([X}) (13) 

where u denotes the growth speed. For homogeneous growth at rate v g in proximal-distal direction y we 
then have u—k-y — — y such that 

f! + \ ([X] + S v,x„ - M + _4_ ( [x| + pffl) . »- m + Rm . 

To model local growth at the limb boundary (Fig. 4) the direction of growth was set to be normal to 
the boundary and proportional to the local FGF concentration, i.e. u = v g F m n where v g is the growth 
speed and n is the unit normal vector on the boundary. The exponent m was set to four, even though 
we note that similar results can be obtained with m = 2. 



Patterning on the limb domain 

Since the model is an example of a Schnakenberg Turing-type model [90] we expected to find parameter 
ranges for which we would observe the emergence of patterns. To judge whether such mechanism could 
yield realistic patterns in a realistic time frame we solved the model on a domain that had been extracted 
from limb bud images taken at Ell. 5. Sox9 expression is induced in response to BMP signaling [91], and 
we were therefore interested in the pattern of the BMP-receptor complex (C = R 2 B) which will mark the 
regions of active BMP signaling. This pattern is necessarily distinct from the pattern of BMP expression 
which remains uniform. In spite of uniform expression of BMP we observe the emergence of pattern in 
the BMP-receptor signaling domains that resemble the Sox9 pattern observed in experiments (Fig 1J-L). 
Thus initially (at r = 150) we observe strong BMP signaling (BMP-receptor binding) in the stalk (Fig 1J) 
(marking the stylopod (s) in the limb bud images in Fig IB or the humerus (h) in Fig IF), two distinct 
spots in the transition zone from stalk to handplate (marking the zeugopod (z) in the limb bud images 
in Fig IB or radia (r) and ulna (u) in Fig IF) and a homogenous expression in the handplate (marking 
the autopod (a) in the limb bud images in Fig IB; digit 4 can be discerned already in Fig IF). Within 
At = 200 the pattern in the handplate splits in the simulation (Fig IK) and distinct condensations 
can be discerned that mark digits 2-5 in the limb bud images (Fig 1C,G). For the horseshoe pattern 
to split the BMP expression rate at the proximal end of the domain (rate pvi) must be 3- fold higher 
than the BMP expression rate in the "hand plate" (rate Pbi)- Without such a local increase in BMP 
expression the horseshoe pattern would not split. With time the digit condensations further elongate 
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(Fig 1D,E,H). In the simulations the digit condensations split while elongating (Fig 1L). Elongation of 
patterns becomes more realistic if we solve the model on a uniformly growing domain (Fig. 3E) or if 
we allow the domain to expand in direction of the highest FGF concentration (Fig 4). Such deforming 
growth towards higher FGF concentrations is also observed in cultured limb buds [38,92]. To achieve 
the characteristic digit shape it is important that the FGF concentrates distally from the BMP-receptor 
spots (Fig 4A-C) as also observed in limb buds (Fig 4D). This is achieved only if BMP-receptor signaling 
has pre-dominantly a positive impact on FGF expression, i.e. if the BMP concentration is relatively low 
such that «i ~ R 2 B -C «2- We therefore predict that at the time of digit formation BMP signaling 
enhances rather than represses FGF8 expression. We should stress that this is necessary only in order 
to obtain the correct positioning of the FGF8 expression domain. The BMP-receptor pattern described 
above can be obtained independently of whether BMP regulates FGF signaling positively or negatively. 

In our simulations 5 spots emerge simultaneously in the handplate while in the developing limb digits 
appear in a sequence (4,2,3/5,1) with digit 1 appearing much later than the other digits. Also the 
appearance of digit 1 appears to be regulated in a different way from the other digits. Thus digit 1 
is still formed in SHH knock-out mutants [10] while digit 1 is lost in FGF8 knock-out mutants [38]. 
The patterns depend on the specific choice of parameters (Table 1), and four instead of five spots could 
easily be obtained by altering any but the receptor expression rate (Fig. 5A). Alternatively the rate of 
FGF expression pp could be lowered 10-fold on the anterior side initially and then increase over time; 
such asymmetry in expression pattern is indeed observed in the limb bud. Thus while Fgf8 is expressed 
uniformly in the AER, expression of Fgf4 appears to be biased to the posterior site of the limb bud. 
FGF8 is expressed already early during limb bud development and remains expressed throughout the 
patterning process while Fgf4 expression is first detected around E10.5 and ceases around E12. As a 
result we expect enhanced FGF production in the posterior part of the limb bud at the onset of Sox9 
expression. However, in this paper we will continue to analyse homogenous expression patterns. 

Regulatory impact of the parameters 

The patterns depend on the specific choice of parameters (Table 1) and since the parameters are difficult 
to determine accurately in experiments we wondered how sensitive our results would be to variations 
in parameter values. We find that all model parameters affect the pattern (Fig. 5A). Nonetheless the 
stable parameter ranges are sufficiently wide that molecular noise should not affect the patterning pro- 
cess. Parameters with particular wide stable ranges include 8b, the rate of BMP degradation, which can 
be varied 6- fold without affecting the pattern, as well as the expression rates pt2, Pr and pp for BMP, 
receptor, and FGF and the threshold for BMP-dependent FGF expression, K\. 

Since we know the time scale of the process we can relate the non-dimensional kinetic rates to their 
dimensional counterparts. In our simulations splitting of the pattern in the handplate takes about 
At = 100. The processes of one embryonic day (86400 seconds) therefore corresponds to At > 200 
such that t = 1 corresponds to at least t = T x r = 432 seconds in dimensional time. Accordingly 
the rate of FGF removal would correspond to <1f = 1/T ~ 2 x 10~ 3 s _1 in dimensional terms which is 
similar to the experimentally reported value of 10~ 3 [87]. Since in our simulations any change in 
the rate of FGF removal can be counterbalanced by a change in FGF expression (Fig. 5B) a two-fold 
lower rate in FGF removal could also be obtained if we reduced the FGF expression rate by about 10-fold. 

The size of the limb bud changes over time but is about 2 millimeters (mm) at E12.5, the final time 
of our simulation. Accordingly the length scale is A = 2mm/ (H + Ri + R2)— 2mm/ 19.35 = 0.1 mm 
and the FGF diffusion coefficient in the model corresponds to 34 pm 2 s _1 . We further require that BMP 
diffuses about 2.7-fold faster than FGF which brings the BMP diffusion constant close to the diffusion 
constant D ~ 100 pm 2 that is typically measured for soluble proteins [87] unless diffusion is impeded 
by adsorption to the surface [86]. We further require that the BMP receptor diffuses about 100- fold more 
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slowly than the BMP protein. The typical range for diffusional constants of membrane proteins is several 
orders of magnitude lower than for proteins in solution (i.e. D = 0.1 — 0.001 /ito 2 s~ 1 [70,77]). A 100-fold 
reduced diffusion constant may still overestimate the receptor diffusion constant, in particular because 
the receptor would then diffuse a distance of 3.5 cell lengths (I — ^/2Z?flt 1 / 2 — 35 /Ltm) within its receptor 
half- life time of t 1 / 2 = ln(2)/dri = 960 s even though receptors are mainly restricted to the surface of 
single cells. We therefore also tested alternative parameterizations in which the receptor diffusion con- 
stant is 400-fold lower than the BMP diffusion constant. Importantly we still obtain qualitative similar, 
yet smaller and sharper, patterns also with a much lower receptor diffusion constant as long as the FGF 
and BMP expression are adjusted accordingly. 

The pattern is also very sensitive to some of the protein expression rates. While there are no reported 
measurements of protein expression rates and concentrations for the limb bud that these constraints could 
be related to, there is a large number of mutant phenotypes that can be used for a qualitative comparison. 
These will be discussed in the following. 

FGF 

Many different FGFs are expressed during limb development. In the model we consider only those FGFs 
that are important during digit formation, i.e. FGF4 and FGF8. In mice that are lacking both Fgf4 
and Fgf8 function in the forelimb AER limb bud mesenchyme fails to survive and the limb does not de- 
velop [ ]. FGF4 is expressed only transiently in the limb bud (between E10 and E12 [37,43]) and removal 
of only FGF4 does not affect digit patterning [39,94]. The effect of Fgf8 removal depends on the develop- 
mental stage. Inactivating Fgf8 in early limb ectoderm (Msx2-cre;Fgf8 hindlimbs and in RAR-cre;Fgf8 
forelimbs) causes a substantial reduction in limb-bud size, a delay in SHH expression, misregulation of 
Fgf4 expression, and loss of digit 1 (the thumb) [38]. In contrast, when functional Fgf8 is transiently 
expressed, as in Msx2-cre;Fgf8 forelimbs, digit 1 is present and digit 2 or 3 is lost [37]. If we remove 50% 
of Fgf expression in our model we observe less separation between the spots, similar to what would be 
expected in syndactyly (Fig. 6A). To lose a digit we need to reduce FGF expression to less than 10% of 
its normal levels. In this case we lose the middle digit as is the case in Msx2-cre;Fgf8 forelimbs without 
causing syndactyly. One explanation for the quantitative differences could be the dominance of FGF8 
expression at later stages of limb bud development such that removal of Fgf8 corresponds to removing 
90% of FGF activity. 

All of the skeletal defects caused by inactivation of Fgf8 can be rescued when Fgf4 is expressed 
in place of Fgf8 [95]. On the protein level FGF4 can thus functionally replace FGF8 in limb skeletal 
development. An increase in FGF signaling that occurs when the Fgf4 gain-of- function allele is activated 
in a wild- type limb bud causes formation of a supernumerary posterior digit (postaxial Polydactyly), 
as well as cutaneous syndactyly between all the digits [95]. Increasing Fgf expression by 50% indeed 
results in supernumerary digits (Fig. 6B). If the additional expression is restricted to the posterior side 
(where Fgf4 is predominantly expressed) then the 50% increase results in the posterior placement of 
one supernumerary digit as well as a merging of spots on the anterior side which would correspond to 
syndactyly (Fig. 6C). 

The model further predicts that the expression and distribution of FGFs remains homogenous initially, 
but subsequently forms distinct patches as indeed observed in experiments (Fig. S3G-J). 

BMP 

There are several BMPs expressed in the developing limb, including most importantly BMP2, BMP4, 
and BMP7 [65]. Much as observed for BMP2 in experiments [79] we predict BMP to be first expressed 
rather homogenously at the boundary of the domain and to be subsequently expressed most strongly 
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outside the digit domains. In our simulation a lower BMP production rate leads to digit loss (Fig. 6E) as 
indeed observed in the BMP2:BMP4 conditional knock-out [66]. A 30% increased production rate results 
in Polydactyly in our simulation (Fig. 6F) and ectopically expressed BMP2 indeed induces duplication 
of digit 2 and bifurcation of digit 3 [ ]. Moreover, conditional inactivation of the gene for BMP4 results 
in Polydactyly [ ] and other experiments reveal that even though conditional inactivation of a gene for 
BMP4 downregulates BMP4 expression, the overall amount of BMP signaling increases [68]. As the BMP 
expression rate is further increased in our simulations (> 50% increase) the patterns merge to give rise 
to polysyndactily (Fig. 6G). Such strong over-expression has not yet been achieved in mutants so that 
we cannot check this prediction. 

BMP receptor 

Conditional mutants of BMP receptor type IA lack digits. Since the receptor is an integral part of the 
proposed Turing mechanism we obtain loss of digits also in the model. Hypomorphic mutants of the BMP 
receptor have not been described; the model would predict loss of digits (Fig. S3). Mild overexpression of 
the receptor should lead to a gain of digits while stronger overexpression will lead to a loss of patterning 
(Fig. S3). 

We conclude that the expression rates of BMP and BMP receptor need to be tightly controlled to 
permit patterning. While FGF is not an essential part of the patterning module FGF enhances BMP 
production and can thus affect the number of digits. 

The impact of domain size and shape 

Turing patterns depend on the size of the computational domain and more modes (patterns) emerge as 
the size of the domain increases [90]. Several mouse mutants have been reported with different shapes 
and sizes of their limb domain, and larger domains typically correlate with increased digit numbers. 
Thus GLI3 mutants are polydactylous and the limb domain seems to be widened to 130% while Gremlin 
mutants are oligodactylous and the limb domain is shrink to about 60% of its original size [10,73]. 
Previous theoretical models showed that such a correlation can in principle be explained with Turing 
models [22], but the domain shapes were not based on real limb buds and no quantitative comparison 
of domain size and digit number was previously carried out. When we shrink the computational domain 
uniformly to 60% of its original size we observe a decrease to two BMP-receptor (BR 2 ) patches (Fig. 
7A,B), much as is observed in experiments [73]. When we enlarge the computational domain uniformly 
to 130% of its original size the digit number more than doubles (Fig. 7A). However, in the GLI3 _//_ 
mutant the expansion is observed only along the anterior-posterior axis [10]. If we increase only the 
anterior-posterior axis of the limb bud to 130% we observe four additional digits which emerge mainly 
in the middle of the domain (Fig. 7C,D). A closer analysis of GLI3 _//_ mutants reveals that limb buds 
are expanded mainly to the anterior side and that extra digits are concentrated in the anterior domain, 
while the posterior side is about normal. When we analysed a domain that was expanded by 160% to the 
anterior and constant on the posterior side we observed three additional spots which now emerge in the 
anterior side (Fig. 7E,F). This corresponds well to the 8 digits observed in the GLI3-/- null mouse [10]. 

The impact of growth 

The limb bud grows while patterns emerge. Given the importance of domain size and shape we wondered 
how growth would affect the patterning process. We measured images of limb buds between E9 and 
E13 and noticed an almost linear growth in proximal-distal direction (Fig. 3A). Limb buds increase in 
length by about 5- fold within the first four days of development. When we solved the model on a linearly, 
uniformly growing domain we obtain similar patterns as before (Fig. 3C). However, the number of digits 
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strongly depends on the growth speed (Fig. 3B). For low growth rates we can adjust the number of 
digits with the help of the other parameters. For increasing growth speeds the speed greatly affects the 
patterning process and we both lose and gain digits (Fig. 3B). As the growth speed further increases 
patterns merge and eventually all digits are lost (Fig. 3B). From Figure 3A we can estimate the growth 
speed as about 7.7 nanometers per second. In dimensionless terms (T = 432 sec and A = 0.1 mm) this 
corresponds to a growth speed of about v g = 0.03 which is too large to permit patterning (Fig. 3B). This 
may be the reason why there is a the complex pre-patterning process that we do not consider in this model 
and which starts already around E9.5 [1,2]. The complex regulatory interactions prior to the emergence 
of digits adjust the expression domains of BMP and FGFs in the limb bud. Such embedding of a Turing 
mechanism into a sophisticated regulatory network may be important to permit the patterning process 
to proceed on a rapidly growing domain. Further experiments combined with more detailed modelling 
also of the other parts of the regulatory network present in the developing limb will, however, be required 
to either confirm a receptor-ligand based Turing mechanism (embedded in a wider network) or to define 
an alternative mechanism for digit patterning on the rapidly growing limb bud. 

Discussion 

Limb development has been studied for decades and much is known about the molecular networks that 
regulate limb development [1,2]. Based on the available information BMP signaling appears to be the 
key signaling pathway that regulates the expression of Sox9, the first marker of digit condensations. 
When modelling BMP signaling we noticed that the regulatory interactions with its receptor constitute 
a Turing mechanisms. Pattern thus easily emerge. The emerging pattern sensitively depended on the 
shape and size of the domain, and only when we solved the set of equations on a domain that had been 
extracted from limb bud images we could reproduce patterns of BMP-receptor activity that resembled 
the experimentally observed patterns for Sox9. An important test for the suitability of a mathematical 
model is its consistency of model predictions with a wide range of independent experimental observations. 
Limb bud development has been studied intensively and a large body of experimental results exists to 
test the model with. If parameters are chosen appropriately also the phenotype of available mutants can 
be reproduced. The mechanism is robust to small changes in the parameters as may arise from molecular 
noise. 

Interestingly, while we achieve good correspondence to available experimental data on the static 
domain, formation of a Turing pattern fails on a domain that grows as fast as measured in the embryo. 
We propose that the early pre-patterning that involves also SHH/Ptch/Gli signaling as well as the many 
other regulatory interactions that we ignore in this simple model are in place to support the patterning 
mechanism on a rapidly growing domain. We are currently developing a simulation for these processes 
to better understand the impact of these earlier regulatory processes. 

Many models have previously been proposed to explain the emergence of digits, and Turing models 
in particular have been found previously to be sufficiently flexible to recapitulate the various patterning 
aspects of digit emergence [30,31,34]. So far, however, it had been difficult to link the components of a 
potential Turing mechanism to the molecular constituents. We propose that the Turing pattern in the 
limb arises from a ligand-receptor interaction. For this to work the ligand needs to be a dimer as is the 
case for BMPs [78]. Moreover, the receptor needs to diffuse, yet at a much lower speed than the ligands 
as is indeed the case. We notice that the diffusion constant that we are using in this simulation is on the 
high end of what is plausible; lower receptor diffusion constants can be accommodated but the pattern is 
then sharper and smaller. Again it is possible and plausible that the additional regulatory interactions 
in the limb bud help to expand the pattern. 

We have previously proposed a mechanism for branch point and branch mode selection in the develop- 
ing lung [96]. Interestingly, also here receptor-ligand interactions gave rise to a Turing pattern that could 
recapitulate the different observed branching patterns in the developing lung for physiological parameters 
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and a similar mechanism can be identified also for the developing kidney (Menshykau et al, submitted). 
While all three systems are based on different signaling proteins (i.e. BMPs in the limb, SHH in the 
lung, and GDNF in the ureter) they may all exploit the same mechanism to pattern an embryonic field 
in that these proteins are all multimers that interact with their receptors in a way that induces receptor 
expression. Interestingly, while BMP can affect its own expression either in a positive or negative feed- 
back via FGFs [52,53], SHH in the lung clearly engages in a negative feedback via FGF10, while GDNF 
in the ureter engages in a positive feedback via Wnt signaling [97]. We can obtain the same pattern with 
a positive and a negative feedback in the limb, lung, and ureter. However, we notice that the possible 
parameter space is wider in case of a positive feedback. It will be interesting to model further embryonic 
self-organizing systems to see whether a receptor-ligand interaction may be a general paradigm to enable 
Turing pattern to emerge in developmental systems. 

Further advancements in our understanding of limb development will require the development of three- 
dimensional models and the inclusion of more signaling factors. The parameterization and validation of 
such models will require new experimental data that also reveal the three dimensional dynamics. Such 
information can now be acquired with the help of optical projection tomographs [ ]. Further advances 
in experimental techniques can thus be expected to provide exciting new insights into the regulatory 
processes of digit formation during limb development. 

Methods 

Numerical Solution of PDEs 

The PDEs were solved with finite element methods as implemented in COMSOL Multiphysics 4.1 and 
4.2. COMSOL Multiphysics is a well-established software package and several studies confirm that 
COMSOL provides accurate solutions to reaction-diffusion equations both on constant [99] and growing 
two-dimensional domains [100-102]. Both mesh and the time step were refined until further refinement 
no longer resulted in noticeable improvements as judged by eye. The simulations were optimised for 
computational efficiency as described in [103]. 
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Figure Legends 
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Figure 1. The emergence of pattern on the limb bud-shaped domain. (A-H) The expression 
pattern of Sox9 in mouse embryos with developmental time (in embryonic days). Panel (A-E) was 
reproduced from Fig 1A-E in [44]; Panel (F-H) was reproduced from Fig 2 in [45]. (J-L) The 
distribution of the BR 2 complex at different simulation times (J) r = 150 , (K) r = 350 ((L) r = 750. 
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Figure 2. The model. A The network represents the regulatory interactions that were considered in 
the model. Thus BMP binds the receptor reversibly to form BMP-receptor complexes. BMP-receptor 
complexes induce the production of receptors and enhances FGF activity. FGF induces BMP 
expression. B The effective regulatory interactions as captured by Eq (9)-(ll) . Thus BMP, B, has a 
positive impact on receptors, R, and on FGF, F, (both via BMP-receptor complexes that are not 
included here). Receptors are auto-activating (when bound by BMP), while BMP are auto- inhibitory 
(as they enhance their own decay by receptor binding). FGF and BMP are mutually enhancing each 
other. For more details see text. (C) The limb bud shape at Ell. 5, adopted from [45]. (D) The 
computational domain. The values for the different lengths are summarized in Table SI. Ri are the 
radial axes of the elliptical limb bud. Hq and Wq are the height and width of the stalk. H\ represents 
the height of the domain in the stalk where the BMP expression is enhanced. The dashed line indicates 
the AER where AER-FGFs are expressed. For details see text. 
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Figure 3. Impact of growth on patterning. (A) The increase in the length of the proximal-distal 
mouse limb bud axis with developmental time (in embryonic days) measured from forty limb buds 
images. Linear growth can be observed at rate of 0.6 ±0.1 mm per day (B) Impact of the number of 
spots (digits) on the growth rate. 5 spots (digits) can be obtained over a wider range of growth rates 
only if we reduce the FGF production rate by 3-fold (pF gr0 wth ~ PF s tatic/3, red line); 
PFgrowth — PFstatici black line). (C-E) Snapshots of simulated limb buds at different points. The 
domain was grown in the proximal-distal axes linearly from 20% of final size in 8000 timesteps with 
PFgrowth ~ PFstatic/3- To prevent the horse shoe patterning from splitting up into digits early, p* B had 
to be decreased from 1.3 to 1 by 0.1 every 2000 time and we set ps2 = 5 • psi to ensure the stalk has 
uniform patterning. (C) Limb bud at r = 4000 and 60% of final size.(D)Limb bud at r = 5300 and 
73% of final size.(E)Limb bud at r = 8000 and 100% of final size. (F-H) The expression pattern of 
Sox9 in mouse embryos with developmental time (in embryonic days) as reproduced from Fig 2 in [45]. 
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Figure 4. Free deformation of the domain according to the local concentration of FGF. 

The FGF patterning gets split up at the positions of the digits due to the positive feedback of the BR 2 
complex on FGF. The domain can then be deformed by allowing growth normal to the surface 
according to the local concentration of FGF. The growth rate is equivalent to [FGF] 4 x v g . (A-C) 

Expression of FGF on the boundary (i.e. ps + p* B (Rvgw^K^pyi ) at different simulation times t = 500, 
1500 and 2500. (D) FGF8 expression pattern in mouse limb at E13.5; adapted from Fig 2C in [95]; 
(E-G) The spatial distribution of the BR 2 complex at r = 500, 1500 and 2500. (H) Sox9 expression 
pattern in mouse limb buds at E13.5; adapted from Fig 8M in [1 n .']. Note that the rates of receptor and 
BMP degradation were increased by 10 percent relative to the standard values used. 
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Figure 5. Parameter dependency of patterning. (A) Local stability analysis. Each parameter as 
perturbed independently from the reference value given in Table 1 as indicated on the vertical axis and 
the range was recorded for which pattern were qualitatively preserved on a static domain. The numbers 
above and below the bar indicate how many digits were either gained (+) or lost (— ) as the range was 
exceeded. (B) The parameter space for which the pattern is preserved is larger (black line) if 
parameters are changed together as illustrated for the FGF production and degradation rates. The 
simulations were analyzed at r = 750. Please note that the FGF degradation rate determines the time 
scale of the simulations (T — 1/dp). Compensation with the FGF production rate thus allow us to 
simulate the processes on different time scales. 
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Figure 6. The impact of protein production rates on digit numbers. (A-D) The impact of 
changes in the FGF production rate on patterning. (A) A 91% reduction in the FGF expression rate 
results in the loss of one pattern (digit). (B) Additional spots emerge as the FGF production rate is 
increased uniformly by 50%. (C) One additional spot emerges on the posterior site (RHS) if FGF 
expression is enhanced by 50% only on the posterior site (RHS). Spots on the anterior side (LHS) merge 
as characteristic for syndactyly. (D) The number of digits at different relative FGF production rates (1 
corresponds to the reference rate in Table 1) for r = 750. Loss of digits at higher FGF production rates 
is due to presence of stripes; spots can be recovered if simulations are run longer. (E-H) The impact of 
changes in the BMP production rate. (E) Digits are lost when the BMP production rate is reduced to 
50%. (F) Additional spots emerge as the BMP production rate is increased by 30%. (G) The pattern 
merge (polysyndactily) as the BMP production rate is increased by 50%. (H) The number of digits at 
different relative BMP production rates (1 corresponds to the reference rate in Table 1) for r = 750. 
Note that at higher production rates pattern merge and digits are no longer observed. The constant but 
spatially modulated BMP production rates psi, PB2 (black and blue lines) and the FGF-dependent 
BMP production rate ps* (red line) have different effects on the patterning. The spatially modulated 
BMP pB2 in the stalk has little impact at larger values because most BMP then binds to receptors in 
the stalk. The spot(s) at the proximal end of the stalk and the two spots appearing at the transition 
from the stalk to the circular domain were preserved under all conditions and thus not counted. 
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Figure 7. Impact of domain size on patterning. (A) Number of spots (digits) on a domain that 
is changed uniformly in its size. (B) Patterning on a domain that is uniformly shrunk to 60% of its 
normal size. (C) Number of spots (digits) on a domain where the length of the posterior and anterior 
axes (R3 and R4) are changed as indicated. (D) Patterning on a domain where the AP axis (R3 and 
R4) is expanded to 130% of its normal size. (E) Number of spots (digits) on a domain where the length 
of the anterior axis (R3) is changed as indicated. (F) Patterning on a domain where the anterior axis 
(R3) is expanded to 160% of its normal size. Note: The spot(s) at the proximal end of the stalk and the 
two spots appearing at the transition from the stalk to the circular domain were preserved under all 
conditions and are thus not counted in panels A,C,E. 
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Tables 



Table 1. Values of the dimensionless parameters. 



Parameter 


Non-dim value 


Description 


Production Rates 




OTP 


243.1 


max production rate of F 


PR 


0.0145 


max production rate of R 


IJ 


1.361 


max C-dependent rate of R production 


< B 




max F-dependent rate of B production 


Pbi 


0.063 


Constant production rate of B in the "handplate" 


PB1 


3 • pbi 


Constant production rate of B in the lower stalk (Hi) 


Diffusion Coefficients 




D B 


2.7 


BMP diffusion constant 


D R 


0.01 • D B 


Receptor diffusion constant 


Decay constants 




Sr 


0.135 


Receptor decay constant 


5b 


0.15 • 5 r 


BMP decay constant 


5c 


u/3 


Receptor-ligand complex decay constant 


Regulation 




n 


2 


Hill coefficient 




0.06588 


Hill constant for BMP -4- FGF positive feedback 




very large 


Hill constant for BMP — FGF negative feedback 


Domain Size 




R Ref 


7.74 


Reference radius 
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Supplementary Text 

Linear Stability Analysis of a simplified version of the proposed model 

We carry out a linear stability analysis to check whether the proposed model is indeed of Turing type as described 
in [1]. The boundary conditions in the proposed model are rather complicated and the spatial modulation of 
boundary conditions and Bmp production introduce further complications (Equation 12 in the main manuscript). 
Here we present the results of a linear stability analysis for a simplified model with constant FGF production on the 
entire domain and no spatial modulation for BMP. Most of the techniques developed for linear stability analysis on 
a growing domain assume exponential growth of the domain (e.g. Madzamuze A 2008 Int J Dyn Sys Dif Eq and 
Madzamudze A et al 2010 Math Biol); similarly we will also assume exponential growth for this linear stability 
analysis. 

The Jacobean J at the steady state in the absence of diffusion can be determined as 



2BR(v - 25 c ) - S r 
-2BR5, 





2BF 2 Rp bl 
(1+F 2 )(1 + BR 2 ) 2 



V 



39.8 172 
-97.7 -194 9.67 
-374 



R 2 S r 



" (1+F 2 )(1+BK 2 ) 2 









2F Pb 



(1+F 2 ) 2 (1 + BR 2 ) (1+F 2 )(1+BR 2 ) 



(1) 



We find that for the parameters used in the manuscript all eigenvalues are negative (i.e. tr ( J) < 0; det( J) < 0) 
in the absence of diffusion, and therefore the steady state is linearly stable. Here we note that the FGF production 
rate that we had to keep constant was varied in a broad range and for the entire range no positive eigenvalues were 
observed in the absence of diffusion. 
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In the presence of diffusion on a constant domain the Jacobean reads as 

( -D r k 2 + 2BR{u - 2S C ) - 5 r R 2 {v - 2S C ) 

-2RRS 2BF 2 R Pbl _n,h 2 -Su- R 2 S F 2 -B 2 Ptll 

z,nno c (i + f 2 )(i+br 2 ) 2 ° b n ° c (i+F 2 )(i+Bfl 2 ) 2 





39.8 - 0.027fc 2 172 

-97.7 -194 - 2.7fc 2 






9.67 

-374 - k 2 



\ 



2F 3 Pti , 2Fp bl 

' (1+F 2 ) 2 (l+Sfl 2 ) (1 + F 2 )(l+Bfi 2 ) 

-D f k 2 - 5 f 



(2) 



J has at least one positive eigenvalue, e.g. k — 20, A = 15.9. Similarly on a growing domain we find A = 15.7 
for k = 20. Therefore we can conclude that the instability is diffusion-driven and therefore of Turing type. 



The linear stability analysis that we carried out here required multiple simplifications and the predicted wavenum- 
ber differs significantly from those observed in the numerical studies. It therefore seems that carrying out a more 
detailed analysis of the simplified system will not bring significant insights for the understanding of the full model. 
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Supplementary Figures 



ABC 




Figure SI. Alternative Implementations of FGF production in the AER. (A) The standard model with FGF 
implemented as flux. (B-E) A thin layer is implemented inside the boundary of the computation domain at a 
distance of 0.02 ■ Ra e f as illustrated in panel B. (C) Thin layer with flux boundary condition. (D) Thin layer 
where FGF production is enhanced by BMP receptor signaling. The simulation was carried out with a slightly 
lower (90%) FGF production rate p s p m — 0.9 pf- (E) Thin layer with constant production of FGF on a thin layer. 
The simulation was carried out with pp 1 " 1 = 0.8pF (F) Eq. 1 1 is defined on the boundary. 
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Figure S2. The spatial distribution of BMP and FGF expression and activity. (A-C) The spatial distribution 
of BMP expression (i.e. ps + P*b ( ) ( r 2 b+i ))' 0^"^) tne concentration of FGF at different simulation times, 
and (G-J) the distribution of FGF expression on the boundary (i.e. pp ^2 > B y^ +K n ^2 B yn +K n ) at different 
simulation times: (left column) r = 150 , (middle column) t = 350, (right column) r = 750. 
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Receptor Production Rate 

Figure S3. The impact of protein production rates on digit numbers Impact of the receptor production rate 
(Pr) on the number of spots appearing in the circular (handplate) domain. The spot(s) at the proximal end of the 
stalk and the two spots appearing at the transition from the stalk to the circular domain were preserved under all 
conditions and thus not counted. 
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Supplementary Tables 



Table SI. GeometryParameters. For an illustration see Figure 2C. RR e f is the reference length scale . 



Parameter 


Value 


Description 


i?3 


RRef 


Length of right semi-major axis of the elliptical limb bud 


i?4 


RRef 


Length of left semi-major axis of elliptical limb bud 


Ri 


RRef 


Length of the semi-minor axis of the elliptical limb bud 


i?2 


-Ri * (1 - WZlARlf* 


Length between limb bud center and stalk 


H a 


0.77* R Ref 


Height of the stalk 


Hi 


0.5* H 


Height of the part of the stalk where BMP production is enhanced (p b2 ^ 0) 


W 


1.41* R Ref 


Width of the stalk 



